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Abstract 

A  rigorous  pseudo  two-dimensional  model  to  simulate  the  cycling  performance  of  a  lithium  ion  cell  is  compared  with  two  simplified  models. 
The  advantage  of  using  simplified  models  is  illustrated  and  their  limitations  are  discussed.  It  is  shown  that  for  1C  or  less  discharge  rates  a 
simple  ordinary  differential  equation  (ODE)  model  can  be  used  to  predict  accurately  the  potential  as  a  function  of  time.  For  rates  higher  than 
1 C,  simplifications  to  the  rigorous  model  are  suggested  that  reduce  the  solution  time  for  the  model. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

The  literature  on  modeling  of  lithium  ion  batteries  is  quite 
extensive  [1-4].  The  first  model  with  two  composite  elec¬ 
trodes  and  a  separator  was  presented  by  Fuller  et  al.  [5] 
Simplified  models  were  presented  by  Doyle  and  Newman 
[6]  to  develop  design  correlations  under  limiting  cases.  Ana¬ 
lytic  expressions  for  the  specific  capacity  against  discharge 
rate  were  also  presented  by  these  authors  [7].  The  model  pre¬ 
sented  by  Fuller  et  al.  [5]  was  extended  by  Ramadass  et  al. 
[17]  to  account  for  the  decay  in  capacity  of  the  cell  with 
cycle  number.  A  side  reaction  leading  to  the  formation  of  a 
film  on  the  surface  of  the  carbon  particles  at  the  anode  was 
proposed  to  occur  during  the  charging  process.  The  poten¬ 
tial  drop  across  the  film  was  expressed  as  a  function  of  the 
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film  thickness,  which  varied  with  time  in  accordance  with 
Faraday’s  law.  The  loss  of  active  material  due  to  the  side 
reaction  and  the  resultant  additional  drop  in  the  anodic  over¬ 
potential  were  used  to  account  for  the  capacity  fade  in  the 
cell.  Further  extensions  to  this  model  were  made  by  Sikha  et 
al.  [8]  that  included  the  change  in  the  porosity  of  the  elec¬ 
trode  material  as  a  function  of  time.  In  all  these  models,  the 
concentration  of  lithium  within  the  solid  phase  was  either  cal¬ 
culated  using  the  superposition  principle  [11]  or  solved  for 
rigorously,  using  a  pseudo  second  dimension  along  the  radius 
of  the  particle.  Since  the  concentration  of  lithium  at  the  par¬ 
ticle  surface  is  the  only  variable  of  interest,  this  methodology 
is  cumbersome  and  time  consuming.  A  very  good  approx¬ 
imation  of  the  concentration  profile  within  the  solid  phase 
was  presented  independently  by  Wang  et  al.  [9]  and  Subra- 
manian  et  al.  [10]  based  on  the  integral  approach  outlined 
by  Ozisik  [11].  In  this  second  approach,  the  concentration 
profile  within  the  solid  particle  is  approximated  by  a  second- 
degree  polynomial  whose  coefficients  are  expressed  in  terms 
of  the  average  concentration  of  lithium  inside  the  particle 
and  the  concentration  at  the  surface.  Thus,  the  need  to  solve 
for  the  concentration  profile  within  the  solid  phase  is  elimi- 
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Nomenclature 

Clj 

specific  area  of  the  porous  electrode 
(m2  m-3) 

Aj 

surface  area  of  the  electrode  (m2) 

brug 

Bruggman’s  coefficient 

c 

concentration  of  lithium  (mol  m-3) 

C$J 

average  concentration  of  lithium  in  the  solid 
phase  of  electrode  (molm-3) 

ci  J 

solid  phase  concentration  of  lithium  at  the 
surface  of  the  sphere  (molm-3) 

DlJ 

diffusion  coefficient  of  lithium  in  the  solid 
phase  inside  electrode  (m2  s-1) 

£>2,eff 

effective  diffusion  coefficient  of  lithium  in  the 
solution  phase  (=D2j£^ug)  (m2  s-1) 

d2 

diffusion  coefficient  of  lithium  in  the  liquid 
phase  (m2  s-1) 

F 

Faraday’s  constant  (C  mol-1) 

fi 

j,  side 

exchange  current  density  for  the  side  reaction 
(Am-2) 

J 

applied  current  (A) 

J 

local  volumetric  current  density  for 
intercalation  reaction  (Am-3) 

Js,j 

side  reaction  current  (A) 

Js,j 

local  volumetric  current  density  for  side 
reaction  (Am-3) 

k 

rate  constant  for  electrochemical  reaction 
(A  m-2  mol  m-3)(1+Q0 

L 

length  of  the  cell  (m) 

Ms  j 

molecular  weight  of  the  side  reaction  product 
(kg  mol-1) 

r 

radial  coordinate  (m) 

Rj 

radius  of  the  particle  (m) 

^SEI 

resistance  of  the  film  (£2  m-2) 

R 

universal  gas  constant  (J  mol  K-1) 

t 

time  (s) 

T 

temperature  (K) 

u e 

local  equilibrium  potential  (V)  (as  described 
in  Appendix  A  of  Ref.  [17]) 

V 

cell  voltage  (V) 

X 

coordinate  across  the  thickness  of  the  cell  (m) 

y 

scaled  radial  co-ordinate  (=  -jpL)  (m) 

Greek 

a 

transfer  coefficients  of  the  electrochemical 
reaction 

£ 

volume  fraction  of  a  phase 

0 

local  potential  of  a  phase  (V) 

0 

over  potential  driving  a  reaction  (V) 

K 

conductivity  of  the  electrolyte  (Sm-1)  (as 
described  in  Appendix  A  of  Ref.  [17]) 

<jeff  effective  conductivity  of  an  electrode  (=  asj) 
(S  m-1) 

a  conductivity  of  the  electrode  (S  m-1) 

pSj  density  of  the  side  reaction  product  (kg  m-3) 

8f  film  thickness  (m) 

Subscript 

1  solid  phase 

2  liquid  phase 

j  =  n  or  p 

p  positive  electrode 

n  negative  electrode 

s  side  reaction  property 

sep  separator 


nated.  While  most  of  these  models  adopt  the  porous  electrode 
approach,  Haran  et  al.  [12]  presented  a  simpler  representa¬ 
tion  of  the  electrode.  This  was  first  presented  for  the  metal 
hydride  system  and  later  extended  to  the  lithium  ion  system 
[13,14].  In  this  model,  each  electrode  is  represented  by  a  sin¬ 
gle  spherical  particle.  This  approach  is  popularly  referred  to 
as  the  single  particle  model. 

Each  of  the  above  approaches  has  some  advantages  com¬ 
pared  to  the  others  as  well  as  demerits.  For  example,  while 
the  porous  electrode  model  has  the  advantage  of  providing 
a  sophisticated  account  of  the  various  physical  processes 
occurring  in  a  battery,  solving  the  model  is  very  time  con¬ 
suming.  The  simplified  single  particle  model  is  orders  of 
magnitude  faster  as  shown  in  this  work,  however,  it  does 
not  account  for  all  the  physical  processes,  for  example,  the 
solution  phase  diffusion  limitations  are  ignored,  thus  the 
validity  of  the  model  is  limited.  Incorporation  of  the  poly¬ 
nomial  approximations  [9-11]  of  the  solid  phase  concen¬ 
tration  of  lithium  in  the  porous  electrode  models  preserves 
their  sophistication  while  simultaneously  reducing  the  solu¬ 
tion  time  [16].  In  simulating  the  performance  of  a  lithium 
ion  cell  over  several  cycles,  the  battery  model  has  to  be 
solved  repeatedly  during  each  cycle.  Ideally,  one  would 
expect  for  such  occasions  that  the  model  is  not  time  con¬ 
suming  and  that  it  will  provide  a  realistic  portrait  of  the 
physical  processes  that  occur  during  cycling.  In  this  context, 
two  simplified  models  are  compared  to  the  rigorous  pseudo 
two-dimensional  model  in  terms  of  their  accuracy  in  simu¬ 
lating  the  cycling  performance  and  the  time  required  for  their 
solution. 


2.  Model  equations 

The  three  approaches  considered  in  this  study  are  the  rig¬ 
orous  pseudo  two-dimensional  model  (or  the  P2D  model),  the 
single  particle  model  (or  the  SP  model)  and  the  porous  elec¬ 
trode  model  with  the  polynomial  approximation  (or  the  PP 
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Fig.  1.  Simplification  of  the  rigorous  two-dimensional  model  (a)  to  the  pseudo  two-dimensional  (P2D)  model  (b)  and  the  single  particle  (SP)  model  (c). 


model).  A  detailed  description  of  the  underlying  assumptions 
in  each  model  and  the  set  of  equations  used  are  elaborated  in 
this  section. 

2.7.  The  pseudo  2D  model 

In  this  approach,  the  equations  developed  by  Fuller  et  al. 
[5]  were  modified  to  include  capacity  fade.  The  discussion 
presented  here  follows  the  model  developed  by  Ramadass  et 
al.  [17].  The  solid  phase  is  assumed  to  comprise  of  identical 
spherical  particles  of  a  predetermined  size  and  diffusion  in 
the  radial  direction  is  assumed  to  be  the  predominant  mode 
of  transport.  The  solution  phase  concentration  and  the  poten¬ 
tials  were  assumed  to  vary  only  in  the  A’  coordinate,  as  shown 
in  Fig.  1(a).  The  reduction  of  the  solvent  (ethylene  carbon¬ 
ate)  is  assumed  to  occur  in  addition  to  the  intercalation  of 


lithium  in  the  carbon  electrode  during  the  charging  process. 
The  entire  contribution  to  the  capacity  fade  is  assumed  to 
be  from  this  reaction.  This  reaction  is  assumed  to  have  a 
constant  open  circuit  potential  of  0.4  V  versus  Li/Li+  and 
the  products  of  this  reaction  are  assumed  to  form  a  film  of 
known  conductivity  over  the  carbon  particles.  As  a  result 
of  this  side  reaction,  a  part  of  the  flux  that  enters  the  car¬ 
bon  particle  is  lost  irreversibly,  resulting  in  a  decrease  of  the 
actual  capacity  of  the  cell.  Further  assumptions  regarding 
the  mechanism  of  capacity  fade  and  their  validity  are  dis¬ 
cussed  by  Ramadass  [20].  The  governing  equations  for  this 
model  are  very  similar  to  those  presented  by  Fuller  et  al.  [5] 
except  for  that,  the  flux  at  the  anode  is  split  into  two  com¬ 
ponents:  one  for  the  intercalation  reaction  (/n)  and  another 
for  the  side  reaction  (/s,n)-  The  thickness  of  the  film  formed 
on  the  carbon  particles  is  related  to  the  side  reaction  flux  by 
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Faraday’s  law  as: 


_  7s?nAfs?n 

dt  <JnPs,nF 

where  8  is  the  thickness  of  the  film.  The  resistance  of  this 
film  results  in  an  additional  contribution  to  the  overpotential 
at  the  anode.  Following  Ramadass  et  al.  [17]  no  side  reaction 
at  the  cathode  is  considered.  Hence,  the  overpotential  (rjn)  for 
the  anode  (the  negative  electrode)  is  given  by: 

/  ,  jjO  (^n  +  Js,n)8 

On  =  01  -  02  -  Un -  (2) 

anK 

while  that  for  the  cathode  is  given  by: 

iiP  =  h-$2-  t/p  (3) 

The  diffusion  inside  the  solid  phase  is  represented  by  the 
Fick’s  laws  as: 


dcUj 

dt 


where  c\j  denotes  the  concentration  of  lithium  inside  the 
solid  in  the  electrode  6f  and  D\j  the  corresponding  diffusion 
coefficient. 

The  solution  of  the  set  of  equations  for  this  rigorous  model 
involves  two  different  length  scales:  the  thickness  of  the  cell 
(L)  is  several  orders  of  magnitude  higher  than  the  average 
radius  ( Rj )  of  the  particles  that  constitute  the  electrode.  Quite 
often  the  numerical  instabilities  that  arise  due  to  the  presence 
of  two  or  more  length  scales  are  circumvented  by  the  use 
of  appropriate  scaling  [15,18,19].  In  the  present  case,  the 
methodology  outlined  in  Appendix  E  of  Ramadass  [20]  was 
employed  and  the  radial  coordinate  was  scaled  as  follows: 

y  =  (5) 

Ki 

and  the  scaled  radial  coordinate  is  termed  as  ‘y’  (See  Fig.  lb). 
In  the  solution  of  this  rigorous  model  (see  equations  in 
Table  1),  the  concentration  of  lithium  inside  the  solid  phase 
(c\j)  is  solved  for  at  each  node  point  along  the  ‘y’  coordinate 
and  the  value  at  y  =  L  is  used  in  the  kinetic  expressions  for 
representing  the  surface  concentration  c\j.  Further  details 
about  this  model  are  presented  in  Appendices  A  and  B  of 
Ramadass  et  al.  [17]. 


2.2.  The  single  particle  model 

In  the  second  approach,  the  single  particle  model  [12-14] 
is  considered.  Each  electrode  is  represented  by  a  single  spher¬ 
ical  particle  whose  area  is  equivalent  to  that  of  the  active 
area  of  the  solid  phase  in  the  porous  electrode.  A  schematic 
of  this  model  is  provided  in  Fig.  lc.  This  model  assumes 
that  the  limitations  posed  by  the  solution  phase  of  the  cell 
are  negligible.  Hence,  the  solution  phase  is  not  considered 
while  developing  the  model  equations.  For  example,  02  are 
set  to  zero  in  Eqs.  (2)  and  (3).  The  assumptions  regarding  the 
side  reaction  are  maintained  in  this  approach  as  well.  This 
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model  is  further  simplified  when  the  concentration  within 
the  sphere  is  approximated  by  a  parabolic  profile  [9-1 1].  The 
solid  phase  concentration  is  represented  by  a  second  order 
polynomial  whose  coefficients  are  expressed  in  terms  of  the 
average  concentration  csj  and  the  concentration  at  the  surface 
of  the  sphere,  c\j.  This  reduces  Eq.  (4)  to  a  first  order  ODE 
and  an  algebraic  equation,  eliminating  the  dependence  on  the 
spatial  variable  V.  The  resultant  equations  are:  [9,10,14] 


dc  2 


dx 


=  0 


At  x  =  Ln  and  Ls 


dx 


dx 


At  x  =  Ln: 


The  initial  condition  to  Eq.  (6)  is: 


At  v  =  Ls : 


cs,j\t=o  —  ci,j,  o  (8) 

where  c\j$  denotes  the  concentration  of  lithium  inside  elec¬ 
trode  at  the  beginning  of  charge  or  discharge.  The  surface 
concentration  c\j  is  then  used  in  the  place  of  cp  j\r=Rj  in  the 
flux  expressions.  For  an  elaborate  derivation  of  Eqs.  (6)— (8) 
see  Subramanian  et  al.  [10]. 


#1.P  Q 

dx 


At  x  =  L: 


2.3.  The  porous  electrode  model  with  the  polynomial 
approximation 


#i,p 

dx 


^app 

°p,eff 


(10) 

(11) 

(12) 

(13) 

(14) 

(15) 

(16) 


The  third  approach  considered  here  involves  the  incorpo¬ 
ration  of  the  parabolic  approximation  [9-11]  into  the  P2D 
model.  The  resultant  equations  of  this  approach  (the  PP 
approach)  retain  the  complexity  of  the  previous  models  based 
on  the  porous  electrode  theory,  but  are  mathematically  sim¬ 
pler.  All  assumptions  presented  for  the  P2D  model  hold  good 
for  the  PP  model  as  well. 

In  addition  to  those,  in  the  PP  approach  it  is  assumed  that 
the  concentration  within  each  spherical  particle  of  each  elec¬ 
trode  can  be  approximated  with  a  parabolic  profile  (similar  to 
the  SP  model)  and  as  a  result,  the  solution  of  Eq.  (4)  at  each 
point  along  the  radial  coordinate  is  now  not  required.  In  the 
set  of  governing  equations  for  the  P2D  model,  the  solid  phase 
concentration  (c\j)  which  was  obtained  from  the  solution  of 
Eq.  (4)  is  now  replaced  with  the  concentration  at  the  surface 
(c\j)  as  given  by  the  solution  of  Eqs.  (6)-(8)  above.  All  other 
equations  for  the  PP  approach  remain  the  same  as  the  P2D 
approach. 

The  complete  set  of  governing  equations  for  each 
approach  is  shown  in  Table  1.  Fig.  1  shows  the  simplification 
of  the  rigorous  two-dimensional  model  to  the  pseudo  two- 
dimensional  model  and  into  the  single  particle  model.  The 
solid  phase  potential  0 i?n  is  arbitrarily  set  to  zero  at  v  =  0. 
The  other  boundary  conditions  obtained  from  jump  material 
and  charge  balances  express  the  continuity  of  flux  across  each 
boundary. 

At  v  =  0  and  x  =  L: 


3.  Simulation  of  the  cycling  performance 

Simulation  of  each  cycle  consists  of  three  steps.  Charg¬ 
ing  is  carried  out  at  constant  current  until  the  cutoff  voltage 
for  charge  (4.2  V)  is  reached  and  then  the  cell  is  charged  in 
the  constant  voltage  mode,  until  the  current  drops  to  50  mA. 
It  is  then  followed  by  discharge  at  a  constant  current  until 
the  cutoff  voltage  for  discharge  (3.0  V)  is  reached.  During 
the  simulation  of  cycling  performance,  the  model  is  solved 
repeatedly  maintaining  all  parameters  at  constant  values  (as 
shown  in  Table  2).  The  fade  in  capacity  with  cycling  is 
brought  about  by  the  additional  side  reaction  flux  (/s?n),  the 
increase  in  film  thickness  over  time  and  the  resultant  change 
in  the  anodic  overpotential.  The  solution  of  the  model  equa¬ 
tions  for  simulating  the  cycling  performance  is  different  from 
that  of  a  single  discharge  curve  in  that  there  arise  numerical 
inconsistencies  between  values  the  variables  solved  for  at  the 
time  between  two  consecutive  steps.  For  example,  the  set  of 
values  for  the  dependent  variables  obtained  at  the  end  of  a  dis¬ 
charge  step  does  not  always  provide  consistent  initial  values 
for  those  variables  at  the  beginning  of  the  next  charge  step. 
This  issue  has  been  discussed  in  detail  by  Wu  and  White  [22]. 
As  recommended  by  these  authors  an  indigenous  initializa¬ 
tion  subroutine  DAEIS  [22]  is  used  to  address  this  numerical 
inconsistency.  The  solution  of  the  model  equations  is  car¬ 
ried  out  using  the  FORTRAN  solver  DASRT  [21].  For  each 
approach,  the  corresponding  set  of  model  equations  presented 
in  Table  1  are  solved  repeatedly  to  predict  the  profile  of  the 
cell  voltage  as  a  function  of  time  for  800  cycles. 
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Table  2 


List  of  parameters  (all  values  are  from  Ref.  [17]  unless  otherwise  stated) 


Symbol 

Anode 

Cathode 

Separator 

o  (S  m-1) 

100 

100 

Sl 

0.49 

0.59 

£ 2 

0.485 

0.365 

0.724 

brug 

4.0 

4.0 

4.0 

D\  (m2  s-1) 

3.9e— 14 

l.Oe— 14 

D2  (m2  s— 1 ) 

7.5e— 10 

7.5e— 10 

7.5e— 10 

k  ((Am-2  molm-3)3/2) 

4.854e— 6 

2.252e— 6 

cfax  (molm-3) 

30555 

51555 

c\,o  (molm-3) 

0.03  x  30555 

0.95  x  51555 

c2,o  (molm-3) 

1000 

1000 

1000 

Rj  (m) 

2e— 6 

2e— 6 

U  (m) 

88e— 6 

80e— 6 

80e— 6 

02,0 

0 

0 

0 

Crcf,s  (V) 

0.40 

Rset  (£2  m-2) 

0.01 

'“side  (Am'2> 

l.Oe— lla 

Msj  (kg  mol  ') 

7.3e— 4 

pSj  (kgm-3) 

2.1e— 3 

k&j  (Sm-1) 

l.Oe— 2a 

Aj  (m-2) 

603.06e— 6 

531. 3e— 6 

a 

0.5 

0.5 

t+ 

0.363 

0.363 

F  (Cmol"1) 

96487 

R  (J  molK-1) 

8.314 

T  (K) 

298.15 

a  Assumed. 


4.  Results  and  discussion 

The  validity  of  the  parabolic  approximation  has  been 
established  by  Wang  et  al.  [9]  and  Subramanian  et  al.  [10], 
and  hence  is  not  repeated  here.  As  shown  in  Ref.  [10],  the 
approximation  is  very  efficient  for  low  to  moderate  rates  of 
discharge.  There  are  two  possible  limitations  that  may  arise 
to  this  approximation.  Under  very  high  rates  of  discharge,  the 
approximation  of  the  concentration  profile  inside  a  spherical 
particle  with  a  second-degree  polynomial  is  not  sufficiently 
accurate.  To  address  this  limitation,  Subramanian  et  al.  [10] 
recommended  higher  order  approximations  of  the  concentra¬ 
tion  profile.  The  second  limitation  arises  from  the  inability  of 
the  approximation  to  capture  the  concentration  profile  at  very 
short  times  accurately.  Wang  and  Srinivasan  [23]  suggested 
an  empirical  correction  factor  to  the  average  concentration 
(cSj)  that  is  obtained  from  Eqs.  (6)— (8)  to  simulate  the  tran¬ 
sient  profile  of  the  concentration  more  accurately. 

Fig.  2a  shows  the  discharge  profile  obtained  from  all  the 
three  models  at  0.2C  as  well  as  1C  rates  of  discharge.  As 
observed  at  rates  as  low  as  0.2C,  both  the  approximate  mod¬ 
els  (SP  and  PP)  provide  a  very  accurate  profile.  This  figure 
validates  the  PP  as  well  as  the  SP  model  independently.  Both 
these  models  compare  well  with  the  rigorous  P2D  model. 
There  is  almost  negligible  deviation  of  the  predictions  from 
the  SP  model  from  those  of  the  porous  electrode-based  mod¬ 
els.  The  absolute  values  of  the  percentage  of  error  between 
the  PP  and  the  P2D  models  as  well  as  those  between  the 
SP  and  the  P2D  models  are  shown  in  Fig.  2b  on  a  logarith- 
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Fig.  2.  (a)  Validation  of  the  SP  and  the  PP  models,  (b)  Comparison  of  the 
error  between  the  approximate  models  (PP  and  SP)  and  the  rigorous  model 
(P2D)  at  various  rates  of  discharge. 


mic  scale.  These  show  that  the  SP  model  is  as  valid  as  the 
porous  electrode  models  for  discharge  rates  as  high  as  1 C. 
As  observed  from  this  figure,  the  error  between  the  PP  and 
the  rigorous  P2D  model  is  of  the  order  of  le— 5  for  the  0.2C 
rate.  This  establishes  the  validity  of  the  polynomial  approx¬ 
imation  even  in  a  sophisticated  porous  electrode  model.  The 
set  of  parameters  used  to  obtain  these  curves  is  presented  in 
Table  2.  Fig.  3  shows  the  discharge  curves  predicted  from  the 
SP  and  the  PP  models  for  fresh  as  well  as  800  cycles  at  the 
0.2C  rate.  The  good  degree  of  fit  between  the  predictions  of 
both  the  models  not  only  reinstates  the  results  from  Fig.  2,  but 
also  justifies  the  claim  that  the  SP  approach  is  quite  success¬ 
ful  in  simulating  the  cycling  performance  of  a  cell  as  much 
as  a  rigorous  model  based  on  the  porous  electrode  theory. 
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Table  3 


Comparison  of  the  PP  and  the  SP  models  with  the  P2D  model 


Rate  of  discharge 

Percentage  error  between  the  PP 
and  the  P2D  models 

Percentage  error  between  the  SP 
and  the  P2D  models 

Cycle  1 

Cycle  800 

Cycle  1 

Cycle  800 

0.2C 

le— 5 

4e— 5 

le— 3 

2.42e— 2 

1C 

0.175 

1.112 

3.404a  3.450b 

6.70a7.12b 

2C 

0.013 

0.135 

59.37a  37.30b 

67.43a  47.05b 

a  With  parameters  from  Table  2. 
b  With  modified  parameters  (see  text). 

The  limitations  of  the  SP  model  as  compared  to  the  porous 
electrode  approach  are  presented  in  Fig.  4.  The  discharge 
curves  at  higher  rates  are  shown  for  cycle  1 .  Despite  the  good 
agreement  with  the  porous  electrode  models  until  the  1 C  rate 
there  is  a  significant  deviation  at  higher  rates.  Table  3  presents 
the  total  error  out  of  the  SP  and  the  PP  approaches  as  com¬ 
pared  to  the  P2D  approach,  for  various  rates  at  cycle  1  and 
800.  The  vast  difference  between  the  profiles  predicted  for  the 
2 C  case  is,  however,  expected.  During  low  rates  of  discharge 
the  limitations  from  the  solution  phase  are  not  as  significant 
as  those  presented  by  the  solid  phase.  Hence,  the  assumption 
that  the  solution  phase  limitations  are  negligible  stated  in 
Section  2.3  is  not  violated  under  these  conditions.  The  solid 
phase  is  treated  alike  in  the  PP  as  well  as  the  SP  model;  so, 
the  difference  between  the  two  models  is  virtually  negligi¬ 
ble  at  the  0.2 C  rates.  However,  at  rates  of  discharge  as  high 
as  2 C,  not  only  are  the  restrictions  posed  by  the  transport  in 
the  solution  phase  important  [24],  but  also  the  dependence 
of  the  kinetic  expression  on  the  solution  phase  concentration 
is  significantly  different  for  the  two  models.  Whereas  the  PP 
model  follows  the  change  in  the  solution  phase  concentra¬ 
tion  via  the  material  balance  in  the  solution  phase,  the  SP 


Fig.  3.  Comparison  of  the  prediction  of  the  cycling  performance  from  the 
PP  and  SP  models  at  0.2C  rate. 


model  assumes  the  concentration  to  be  constant  throughout 
discharge.  This  is  reflected  in  the  difference  in  the  slope  of 
the  discharge  plateaus  between  the  two  models  in  Fig.  4.  Note 
that  no  such  difference  is  observed  between  the  two  models 
at  the  0.2 C  rate  in  Fig.  1,  which  substantiates  the  assumption 
that  the  contribution  of  the  solution  phase  is  not  significant 
at  rates  of  discharge  below  1 C.  To  validate  the  argument  that 
the  difference  between  the  PP  and  the  SP  models  at  higher 
rates  is  because  of  the  solution  phase  limitations,  the  solution 
phase  conductivity  (/ceff)  was  increased  by  a  factor  of  two  and 
the  solution  phase  diffusion  coefficient  (D2,eff)  was  set  to  a 
very  high  value  (le— 3  m2  s_1)  in  the  PP  model.  The  revised 
profiles  are  also  shown  in  Fig.  4.  As  observed,  at  1C  rate 
the  change  is  minimal.  However,  at  the  2C  rate,  the  revised 
profile  shows  a  better  agreement  with  the  SP  model;  still  a 
complete  agreement  between  the  two  models  is  not  obtained, 
since  the  value  for  the  solution  phase  concentration  in  the 
kinetics  expression  continues  to  be  a  constant.  The  signif¬ 
icance  of  the  change  in  the  solution  phase  concentration  is 
illustrated  in  Fig.  5.  Whereas  for  discharge  rates  up  to  1C  the 
change  in  the  solution  phase  concentration  is  less  than  10%, 
for  the  2C  rate  it  is  as  high  as  40%.  The  SP  model  completely 


Time  (s) 

Fig.  4.  Comparison  of  the  discharge  curves  predicted  from  the  PP  and  SP 
models  at  higher  rates:  (a)  with  parameters  from  Table  2  (b)  with  modified 
parameters  (see  text). 
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Fig.  5.  Change  in  the  solution  phase  concentration  at  various  rates  of  dis¬ 
charge  as  predicted  by  the  PP  model. 

ignores  the  changes  in  the  solution  phase  concentration,  and 
hence  is  incapable  of  predicting  the  discharge  profile  in  the 
kinetics  dominated  regime.  An  alternate  discussion  is  pro¬ 
vided  by  Botte  and  White  [25],  wherein  the  inclusion  of  the 
activity  coefficient  term  in  the  solid  phase  was  suggested  to 
improve  the  predictions  at  high  rates. 

The  computational  time  for  800  cycles  for  all  the  three 
approaches  is  presented  in  Table  4.  Hundred  node  points 
were  used  to  approximate  each  coordinate  using  the  three- 
point  finite  differences  scheme  for  each  electrode  in  the  P2D 
and  the  PP  models.  The  SP  model  has  a  distinct  advantage 
over  the  porous  electrode  models,  in  terms  of  the  compu¬ 
tational  time,  despite  its  limitations  in  terms  of  accuracy  at 
higher  rates.  Nevertheless,  the  PP  model  exhibits  a  good  effi¬ 
ciency  both  in  terms  of  accuracy  as  well  as  computational 
time,  as  observed  from  Tables  3  and  4.  The  incorporation 
of  the  parabolic  approximation  [9-11]  reduces  the  computa¬ 
tional  time  by  at  least  70%.  The  choice  of  DASRT  coupled 
with  DAEIS  is  found  to  be  a  robust  solver  suitable  for  solv¬ 
ing  the  models  repeatedly  during  the  simulation  of  the  cycling 
performance  of  the  cell,  and  hence  is  ideal  for  capacity  fade 
analysis.  DASRT  [21]  is  a  variable  time  step  solver  with  an  in 
built  zero  crossing  detection  feature.  This  feature  of  DASRT, 
along  with  appropriate  constraints,  enables  one  to  interpo¬ 
late  to  the  exact  time  of  transition  from  one  step  to  another 

Table  4 


Comparison  of  the  computational  cost  for  the  three  models 


Model 

CPU  time  for  800  cycles  (s) 

P2D 

1.12e5 

PP 

9600 

SP 

27 

(charging  to  discharging,  for  example).  The  capabilities  of 
DAEIS  are  addressed  by  Wu  and  White  [22] . 


5.  Conclusion 

Two  approximate  models  (SP  and  PP)  were  used  to  sim¬ 
ulate  capacity  fade  in  lithium  ion  batteries  and  compared  to 
a  complete  model  (P2D).  Both  of  the  approximate  models 
compared  well  with  the  P2D  model  up  to  1C  rate  of  dis¬ 
charge.  Incorporation  of  the  parabolic  approximation  for  the 
solid  phase  concentration  of  the  diffusing  species  signifi¬ 
cantly  reduces  the  computational  time  as  compared  to  the 
P2D  model  for  both  the  SP  and  the  PP  models,  while  still 
retaining  sufficient  accuracy  as  reported  by  Subramanian  et 
al.  [10],  even  for  a  complex  sandwich  model.  The  low  solu¬ 
tion  time  as  well  as  a  high  degree  of  accuracy  identifies  these 
models  as  suitable  for  simulating  the  cycling  performance  of 
a  cell,  when  the  model  equations  are  solved  over  and  over. 
The  SP  models  neglects  solution  phase  limitations,  and  hence 
is  not  accurate  at  rates  beyond  1 C,  when  concentration  gradi¬ 
ents  in  the  liquid  phase  become  the  limiting  factor.  While  the 
SP  model  can  be  used  successfully  to  simulate  the  cycling  of 
lithium  ion  batteries  up  to  the  1 C  rate,  the  PP  model  is  best 
suited  for  rates  higher  than  1 C,  since  there  is  no  significant 
loss  of  accuracy  and  at  the  same  time  the  computational  time 
is  greatly  reduced.  DASRT  is  identified  as  the  suitable  choice 
of  solver  due  to  its  robustness  in  solving  the  model  equations 
over  800  cycles. 
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